A nonlinear evolution equation for sand ripples based on geometry and conservation 
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From geometry and conservation we derive two nonlinear evolution equations for sand ripples. In 
the case of a strong wind leading to a net erosion of the sand bed, ripples obey the Benney equation. 
This leads either to order or disorder depending on whether dispersion is strong or weak. In the 
most frequent case where erosion is counterbalanced by deposition, we derive a new one-parameter 
nonlinear equation. It reveals ripple structures which then undergo a coarsening process at long 
times, a process which then slows down dramatically with the growth of the ripple wavelength. 

PACS numbers: 81.05 Rm 
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Ripples on sand in desert and see are fascinating pat- 
terns jjj]. Despite the fact that sand, and granular me- 
dia in general, are very familiar, in principle, to anyone, 
the understanding of their static and dynamical prop- 
erties still continues to pose a formidable challenge to 
theoretical modeling Q. A continuum description, as 
those well established for simple fluids (Navier-Stokes), or 
solids (Hookes law), is lacking, most likely due to a strong 
interaction of disparate scales. Whatever complex a con- 
tinuum theory might be, it should be compatible with 
symmetries and conservation laws. The present paper 
deals with derivations of two generic nonlinear equations 
for sand ripples based on geometry and conservation. 

A model describing the physical origin of ripple forma- 
tion has been presented in the early forties by Bagnold 
|J . This model is based on energetic saltating grains im- 
pacting the ripple. Since then, several contribution both 
analytical |3|jj] and numerical have allowed further 
elucidation of the problem. A derivation of a continuum 
generic nonlinear evolution equation of sand front is lack- 
ing, however. This means in particular that the question 
of whether ordered, or disordered pattern, would prevail 
is to date unanswered, ft is the main goal of this Let- 
ter to address these questions on the basis of geometry 
and conservation. In the most interesting case where ero- 
sion is counterbalanced by deposition on the average, we 
show here that the generic equation in one dimension 
(when the ripple is translationally invariant in the y— di- 
rection) takes the following form close to the instability 
threshold 
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where h is the ripple profile (in a dimensionless form) 
which is a function of x and t. As seen below this equa- 
tion can always be reduced to one dimensionless parame- 
ter denoted here as v. This equation reveals ripples. The 
wavelength is first dominated by that of the linearly most 



unstable mode. At long time they coarsen before reveal- 
ing a dramatic slowing down of the wavelength increase. 
Note that the equation is local, which is a consequence 
of the proximity of the instability threshold (see below). 

For ease of presentation, we consider a one dimensional 
front. The most natural way of representing a front is to 
use intrinsic coordinates, namely the curvature k(s) as 
a function of the arclength s. Each point of the front is 
labelled by its vector position r(a , t), where a is a time- 
independent parametrization of the curve which can be 
taken at liberty to lie in the interval to 1. We obtain 
for the evolution of the arclength s [lol 
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This is simply obtained by setting ds = ^fgda, where g is 
the induced metric, and integrating s — J yfgda by parts. 
Here we have used the relation dg/dt — 2g(^j- + KV n ), 
where Vt and v n are the tangential and normal velocities 
respectively. In order to derive the evolution equation for 



k, we first need to evaluate the commutator [ 



obtain 
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Applying this identity to 9 (the angle between the vertical 
axis and the normal vector), we arrive at 
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Equations (||) and (||) constitute the evolution equations 
for the arclength and the curvature. These equations are 
general and only geometrical concepts are evoked. 

The tangential velocity is a gauge. Indeed if one con- 
siders the front at some time and its state at later time 
there is no way which allows us to state if a point a (lying 
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on the curve) at time t has becomes point a' at t + At, 
or a" (obtained by displacing a' tangentially) . The tan- 
gential velocit y is fixed once the parametrization of the 
curve is given jlfj]). In contrast, v n is a physical quantity. 

A way of viewing that the tangential velocity disap- 
pears from the evolution equation is to work at given s, 
and not at given a. Using (Q), the curvature evolution 
equation (|j) takes the form 

8k, , d 2 21 OK f s , . 
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Gauges usually introduce nonlocality. Similar formula- 
tions were used in other contexts pj]Jl^| . This equation 
constitutes a general nonlinear equation for curvature dy- 
namics if the normal velocity is known. This formulation 
is powerful in numerical studies. If v n is known as a 
function of geometry then Eq. (||) is a closed equation 
for curvature dynamics. Below we shall write down the 
expression for v n inferred from symmetries and conser- 
vation. We shall then use Eq. @ as a starting point of 
our derivation of evolution equations for the ripples in 
the weakly nonlinear regime in terms of Cartesian coor- 
dinates where the front is represented by z = h{x,t). 

Let us first illustrate our analysis on a situation that 
leads to a well known equation in the literature. This 
will serve later to determine relevant nonlinearities. For 
a rotationally invariant system, the normal velocity must 
necessarily be a function of only those quantities which 
are invariant under any surface reparametrization. In 
one dimension the only quantity which is intrinsic is the 
curvature and its odd derivatives with respect to the ar- 
clength. Generically, the normal velocity has thus the 
form 

v n = C + a\K + a,2K 2 + a^K 3 + bin ss + .... (6) 

From now on most of differentiations will be subscripted 
for brevity. Before proceeding further an important re- 
mark is in order concerning the assumption of locality 
in Eq. (g). For sand ripples discussed below, there are 
two kinds of grains that contribute to the development of 
ripples |^-||] : the saltating ones (traveling on length scale 
l s ~ A, where A is the ripple wavelength) that have high 
kinetic energy, and the the low-energy splashed grains 
traveling in reptation (or hopping) on a scale a which is 
several (typically 6—10) times smaller than A. The saltat- 
ing grains are accelerated by the wind and this provides 
the driving force that governs the motion of the surface 
in aeolian ripple growth and in ripple translation. In the 
most frequent case where erosion is counterbalanced by 
deposition, the population of saltating grains remains al- 
most constant (as recognized already by Bagnold pj and 
in l^-fl). The saltating grains serve merely to bring en- 
ergy into the system, the saltating population exchanges 
almost no grains with the reptating population. The in- 
formation on the surface profile is propagated only by 
reptating grains. Saltating grains loose their memory in 



the course of their flight due to collisions and the tur- 
bulent airflow. It follows then that a (reptation length) 
is the natural candidate |^-|5| for a characteristic length 
scale of interaction. Since a/ A -C 1, the local assump- 
tion is legitimate. The ratio a/ A serves here as the small 
expansion parameter. 

The constant C in Eq. (^) expresses the fact that 
the front advances on the average at constant velocity. 
For a straight front the velocity is C (precisely as what 
happens when a front is collecting particles from outside 
when exposed to a given flux). For a weakly curved front 
we have approximately k ~ ~h xx . Using this in (^3j) and 
upon substitution into (|5|) we obtain 

C 

h t = C - a\h xx - bih xxxx + —hi , (7) 

where we have truncated the expansion to leading order, 
as explained below. First we can have a\ > or a% < 0. 
In the former case the straight front is unstable. This is 
what happens in a large number of situations (see Ref. 
Jl3|| ) . The above derivation has concentrated on the sit- 
uation close enough to the threshold where a\ is small ( 
a\ changes sign upon variation of some control parame- 
ter in a given system). For the truncation to make sense, 
a\ must be of order e (some appropriate small quantity) . 
This is the case sufficiently close to the instability thresh- 
old. For the fourth derivative to be of the order of the 
second one, we must have eh this implies in 

Fourier space that q (the wavenumber) must be of order 
y/e. In other words our equation is expected to be valid 
for slowly varying modulations on the scale of the typ- 
ical length of interest in a given problem. This means 
that modulations in physical units occurs on scales of 
order l/y/e. Now in order that the linear part coun- 
terbalance the nonlinear one, we must have eh xx ~ h x . 
Since x ~ this implies that h ~ e. Using the same 

argument one arrives at the fact that time scales as 1/e 2 
(or equivalently the frequencies of interest are of order 
e 2 ). Once the scalings are known, it is a simple matter 
to show that other permissible nonlinearities (e.g., h xx ) 
are of higher order contribution (i.e., of order e 4 ). 

Note that the constant term in Eq. (Q) is unimportant 
since it can be absorbed on the l.h.s. upon a transforma- 
tion h—>h — Ct. Note also (see below) that b% is generi- 
cally positive in order to ensure a well behaved solution at 
a short scale. The sign of C is however unimportant since 
changing it would simply correspond to the transforma- 
tion h — > —h. The equation can be made free of param- 
eter upon appropriate rescaling. Equation (H) is known 
under the name of Kuramoto-Sivashinsky (KS) . It 

models pattern formation in different contexts (see Ref. 
|]l3f ) . For large enough extent of the front in the x direc- 
tion, it reveals spatio-temporal chaos. If a\ < , there is 
no instability, and there would thus be no need to keep 
the fourth derivative in Eq. ([?]) . Adding a small stochas- 
tic force to C (like shot noise in Molecular Beam Epitaxy 
MBE), we obtain the well known Kardar-Parisi-Zhang 
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|L6| equation, introduced to model kinetic roughening in 
MBE. 

Having illustrated our study on a reference example, 
we are now in a position to deal with ripple formation 
under wind blow. On the one hand the front is not ad- 
vancing on the average, in principle. Thus the constant 
C must be set to zero (this is unimportant as seen below) . 
On the other hand, the wind causes the normal front ve- 
locity to be orientation-dependent. We first concentrate 
on the situation where there is a strong erosion, so that 
the front is surrounded by an atmosphere of flying grains 
(a sort of reservoir). The front motion can thus globally 
loose or gain grains from the atmosphere, so that no con- 
straint must be imposed. In the presence of anisotropy 
(due to the wind) the most natural way is to write 

v n — &iK + + biK ss + ai sin 9 + a^sin 9) 2 

d 2 d 2 
+/9i 5-5 (sin 9) + ft — (sin 2 6) + j lK sin + (8) 

The terms in sin (9 express the fact that the growth ve- 
locity depends on the local slope of the front. In the 
present case, the direction perpendicular to the z axis 
(from which 9 is measured) is favored. Had we wanted 
to give a greater importance to the x direction, we would 
then have expanded v n in power of cos 9. However, it is 
a simple matter to realise that this is unimportant for 
our purposes. In the limit of a weakly curved front, the 
substitution of (||) into (||) yields to leading order 

h t = -a.\h x - aih xx - Pih xxx - b\h xxxx + a 2 h 2 x (9) 

This equation is known in the literature under the name 
of Benney equation |17) . It has been derived recently 
from a microscopic model in the context of step-bunching 
dynamics during sublimation of a vicinal surface jl8],|l9| . 
The same equation arises in other contexts such as phase 
dynamics for traveling modes po| , and in some models 
of traffic flow j2l| . The first derivative term can be ab- 
sorbed in the temporal derivative by means of a Galilean 
transformation (i.e. x — > x — ait). The scaling of space, 
time and amplitude with e are obviously the same as for 
the KS equation. The third derivative (not present in the 
KS equation) is of higher order 1/e 1 / 2 if all the scales in 
the KS part are set to one (as seen before space scales as 
l/^/e, and h ~ e, so that a\h xx ~ h xxxx ~ h 2 ~ e 3 , while 
h xxx ~ e 5 / 2 ). Thus the expansion would make sense in 
principle only if (3\ were small enough (of order y/e), a 
demand whose realization depends on the system under 
consideration. This apparent difficulty can be circum- 
vented by noting that h xxx contributes to the imaginary 
part of the linear growth rate (if we seek a solution in the 
form e lqx+UJt j where ui is the growth rate), which concerns 
thus propagative terms, whereas h xx and h xxxx produce 
real contributions. Thus one can 'split' time into a slow 
part corresponding to growth of perturbations, and a fast 
part corresponding to propagation. 

Upon rescaling of the space, time and amplitude, the 
Benney equation can be rewritten in a form in which only 



one parameter survives. Therefore all the coefficients can 
be set to unity except one, let say j3\. Depending on the 
strength of this coefficient, the dynamics is either chaotic 
for Pi of order or smaller than one (we recover the KS 
dynamics), or exhibits a rather ordered structure drifting 
sideways for /3i of order few unities ]l9| . 

For sand ripples, it seems that in most situations an 
equilibrium between erosion and deposition sets in due to 
the drag force of the transported grains exerted on the 
wind (the greater the number of transported grains is, 
the weaker the wind gets and the less it erodes the sand 
bed). In other words, the number of transported grains 
remains constant in average; there is neither net erosion 
nor deposition. In order to treat that case we should 
impose the global conservation condition. In order to 
ensure this all the homogeneous terms in (^) must be left 
out, except the linear terms in k and sin 9. Indeed the 
area A bounded by the front and some horizontal axis 
behaves in course of time as dA/dt = J dsv n , and as 
global conservation imposes dA/dt = 0, all terms giving 
a non-zero contribution to that area should be removed 
(the sand front moves because either particles have left 
the region of interest, or other grains have landed from 
a neighboring part). In this case, the front dynamics is 
governed to leading order by 

Q 2 2 

ht — —ot\h x — a\h xx — flihxxx — b\h xxxx + ^ 2 ^2 (^u) 

(10) 

For an instability a\ must be positive, b\ must be pos- 
itive as well in order to introduce a short wavelength 
cut-off, while the sign of @2 is unimportant, since it can 
be changed upon the transformation h — ► — h. Obviously 
a.\h x can be absorbed in h t via a Galilean transformation, 
and the sign of (3\ is unimportant as well. Space and time 
scale with e as in the KS and Benney equations, while 
the scale of h here is of order one, in a marked contrast 
with the KS limit. This scaling has an important con- 
sequence, to be discussed below. What makes the other 
higher contribution small is precisely the scaling of space 
(higher and higher derivatives are of smaller and smaller 
contributions). After rescaling (and absorbing h x in ht) 
only one parameter survives and the principal equation 
for ripples can be written in a canonical form (Eq. (0)). 
The linear dispersion relation of Eq. (|]J) takes the form 
(by looking for perturbations in the form h ~ e *qx+wt^ 
uj = q 2 — q A — ivq z . The structureless state is linearly un- 
stable against perturbations with wavenumbers smaller 
than q c — 1. A physical origin of the instability was put 
forward long time ago by Bagnold Q. Here we assume 
that the threshold has been reached, and thus we take a 
negative sign in front of the second derivative. Numeri- 
cal solutions for sizes L > 10A C = 2tt of Eq. (|l|) reveals 
an evolution towards a steady state with a given wave- 
length. In a marked contrast to the KS equation which 
exhibits spatiotemporal chaos (or Benney equation when 
dispersion is small) , the new evolution equation (||) leads 
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FIG. 1. The ripple profile at different time. The consecu- 
tive snapshots have been shifted upward to show the drift. 



to ripple pattern (Fig.l). The wavelength is first close 
to that of the most dangerous mode. At long time, the 
structure coarsens producing thereby wider and wider 
dunes [^2| . Then the coarsening slow down dramatically. 
This feature agrees with experiments An extensive 
study will be presented in the future. 

Finally it is important to make some important com- 
ments, (i) Since h is of order one, while x is of or- 
der l/yfe, the slope is of order y/e, and thus remains 
small (see Fig.l). Had the slope been of order one (as in 
some nonlinear equations pi), avalanches would mani- 
fest themselves, and no steady-solutions would have been 
possible, (ii) Figure 1 corresponds to equation (|]). As 
stated above changing the sign in front of the nonlin- 
ear term corresponds to making an up-down operation 
on Figure 1. Inspecting several examples of sand rip- 
ples conveys the strong impression that it is the situa- 
tion in Figure 1 which seems likely, (iii) Since h ~ 1 and 
A ~ the ratio of the amplitude to wavelength close 

to the instability threshold scales as e 1 / 2 . Thus close to 
threshold the amplitude is several times smaller than the 
wavelength (usually ripples have an amplitude which is 
approximately 10 times smaller than their wavelength), 
(iv) We have deliberately, for sake of simplicity, consid- 
ered a one dimensional structure (that is the ripples are 
translationally invariant along the y direction). The ex- 
tension to two dimensions is feasible, and is currently 
under investigation |^3f . This should be crucial with re- 
gard to the study of possible secondary instabilities of 
ripples. 

In summary we have derived a generic nonlinear equa- 
tion to describe sand ripple dynamics. The study is based 
on geometry and conservation. While apparently close 
equations (e.g., the KS equation) lead to spatiotempo- 
ral chaos, the new equation reveals rather steady ripples 
which coarsen with time. The advantage of this study 
lies in the fact that no matter how complex the physics 
might be, the equations close to the instability thresh- 
old must be of the sort given here as long as symmetries 



and conservations are preserved. In turn, geometry and 
conservation can not, by their very nature, provide the 
values of coefficients The same holds for the existence 
of instability. Thus there is a need in future to derive 
our equation starting from a given 'microscopic' model 
and to determine the dimensionless coefficient v in terms 
of physical quantities. Other remarks are in order. We 
have limited ourselves to leading order, which is valid 
close enough to the instability point p^ . Expansion to 
higher order is straightforward. We have assumed that 
the normal velocity is both local in space and in time. 
This is valid close to the instability threshold. If repta- 
tion length a remains the relevant length scale, we may 
expect our assumption of locality to hold at arbitrary dis- 
tance from threshold, since a/ A is the small parameter of 
the expansion. Finally, it goes without saying that the 
present study can have impact on other systems than the 
sand ripples. 
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